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There  are  a  number  of  algorithms  for  calculating  the  least  squares 
isotonic  regression  function  and  six  of  these  algorithms  are  discussed  in 
Section  2.3  of  Barlow ,  Bartholomew,  Bremner  and  Brunk  (1972).  The  most 
widely  used  of  these  algorithms  is  the  "pool  adjacent  violators  algorithm" 
which  is  applicable  only  in  the  case  of  a  simple  linear  ordering  or  an 
amalgamation  of  simple  linear  orderings.  Cran  (1980,  Algorithm  AS  1U9) 
gives  a  routine  for  computing  the  isotonic  regression  for  the  case  of  a 
simple  linear  ordering  using  the  "Up  and  Down  Blocks"  algorithm.  In  many 
applications  of  isotonic  regression  we  have  more  than  one  independent 
variable  and  the  regression  function  is  restricted  to  be  monotone  in  each 
independent  variable.  Severed  of  the  algorithms  given  in  Barlow  et  ed. 
are  applicable  to  this  restriction  but  those  which  are  reasonable 
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to  program  often  lead  to  excessive  computation  time.  Dykstra 

and  Robertson  (1982)  give  an  algorithm  for  computing  the  least  squares 

regression  function  which  is  increasing  in  each  of  several  variables. 

This  algorithm  uses  successive  one-dimensional  smoothings  and  is  very  effi¬ 
cient,  relative  to  other  available  algorithms.  A  Fortran  IV  program  to 
implement  this  algorithm  for  two  independent  variables  is  given  below. 


STRUCTURE 

SUBROUTINE  SMOOTH  (NROW,  NCOL,  X,  W,  NCYCLE,  G,  ICYCLE,  IF AULT ) 
Formal  Parameters 


NROW 

Integer 

input : 

the  number  of  rows 

NCOL 

Integer 

input : 

the  number  of  columns 

X 

Real  Array  (NROW, NCOL) 

input : 

the  original  values 

W 

Real  Array  (NROW, NCOL) 

input : 

the  original  weights ,  all 

zero 

or  positive 

NCYCLE 

Integer 

input : 

the  maximum  number  of  cycles 

if  the  array  does  not  converge 

GCHECK 

Real  Array  (NROW, NCOL) 

workspace 

WT 

Real  Array  (NROW, NCOL) 

workspace 

G 

Real  Array  (NROW, NCOL) 

output : 

the  smoothed  values  if  con- 

vergence  was  attained,  or  the 
values  after  the  maximum  number 
of  iterations 


ICYCLE  Integer  output:  the  actual  number  of  cycles 

the  smoothing  process  required 
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IFAULT  Integer  output:  a  fault  indicator,  equal  to 

1  if  convergence  was  not  attained 
in  the  specified  number  of 
cycles 

2  if  any  element  of  W  is  nega¬ 
tive 

3  if  rows  and/or  columns  exceed  20 
k  if  a  zero  weight  was  replaced 

by  0.00001 
0  otherwise 

Auxiliary  Algorithm 

SMOOTH  uses  the  subroutine  PA V  (K,X,I0RDER,W,FINALX)  to  apply  the  pool 
adjacent  violators  algorithm  for  the  isotonic  regression  on  rows  and 
columns . 

Constants 

The  constant  EPS  is  used  to  determine  if  convergence  has  been  attained. 

If  each  entry  of  the  array  after  the  current  smoothing  does  not  differ 
from  the  corresponding  entry  from  the  previous  smoothing  by  more  than 
EPS,  the  array  is  considered  to  be  smoothed. 

Restrictions 

The  X  and  W  arrays  cannot  have  more  than  20  rows  and/or  columns. 
However,  if  one  wished  to  smooth  a  larger  array,  the  changes  to  the  sub¬ 
routine  would  be  minimal. 

The  W  array  cannot  contain  any  negative  values.  If  zero  values  are 
encountered,  they  are  replaced  by  10~^  for  the  amalgamation. 
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Time  and  Accuracy 

For  the  9X9  array  in  Table  1+  in  Dykstra  and  Robertson  (1982),  convergence 
was  attained  in  47  iterations  at  a  CPU  time  of  .26  seconds. 

The  constant  EPS  specifies  the  change  from  one  cycle  to  the  next  for  con¬ 
vergence.  It  has  been  given  an  initial  value  of  10  however  if 
greater  or  lesser  accuracy  is  required,  the  value  may  be  changed. 

Related  Algorithms 

Cran  (1980)  gives  an  algorithm  (AS  1U9)  for  amalgamation  of  the  means  in 
the  case  of  simple  ordering,  which  can  be  used  as  an  alternative  to 
SUBROUTINE  PAV. 
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DRIVER  PROGRAM  FOR  SMOOTH  SUBROUTINE. 

DIMENSION  X!20*2C>*  WT!20,20>,  G!20,20> 

READ  !  5*  10  )  NROW*  NCul,  NCYCLE 
0  FORM  AT  !  212*  15  ) 

DO  30  I  =  1*  NR  0 V 

READ !  5,  20  >  C  X!I,J>,  J  =  1*  NCOL  ) 

C  FORMAT  C  10F8.4  ) 

0  CONTINUE 

DO  40  I  =  ]*  NROW 

READ  C  5*  20  )  (  WT(I,J>*  J  =  1  *  NCOL  ) 

C  CONTINUE 

ECHO  CHECK  THE  INPUT  ARRAYS. 

URITEC  6*  50  ) 

0  FORMAT!  1 H 1  ) 

WRITE(  6*60) 

0  FORMAT (  4 X *  34HTHE  ORIGINAL  MATRIX  TO  PE  SMOOTHED,//  ) 

DO  80  I  =  1*  NR OU 

WRITE!  6.  70  )  (  X(I«J),  J  =  1*  NCOL  ) 

0  FORMAT!  4X*  10F12.6  ) 

0  CONTINUE 

WRITE!  6*  90  1 

0  FORMAT!  Ill 4X ,  28HTHE  ASSOCIATED  WEIGHT  MATRIX,  //  ) 

DO  100  1=1,  NROW 

WRITE!  6,  70  >  !  WT(I,J),  d  =  1,  NCCL  > 

00  CONTINUE 

SMOOTH  THE  X  ARRAY. 

CALL  SMOOTH!  NROW,  NCOL,  X,  WT,  NCYCLE,  ICYCLE,  IF  AULT,  G  > 

CHECK  FOR  CAULT  ERRORS  IN  THE  SUBROUTINE. 

IF  !  I F  AULT  .EO.  2  >  GOTO  160 

IF  (  I F AULT  .EQ.  3  >  GOTO  180 

IF  !  IF  AULT  • E3 .  4  )  WRITE!  6,  110  > 

10  FORMAT  <  / /4  X,  3 1HZERO  WEIGHTS  HAVE  BEEN  REPLACED,/, 

*4X,  24HBY  0.00001  FOR  SMOOTHING  ) 

OUTPUT  THE  SMOOTHED  MATRIX. 

WRITE!  6*  5  C  ) 

WRITE!  6,  120  )  ICYCLE 

120  FORMAT!  // 4 X,  35HTHE  NUMBER  OF  CYCLES  COMPLETED  WAS*  15  > 

IF  !  I F  AULT  .EO.  1  >  WRITE!  6,  130  > 

130  FORMAT!  4 X ,  4  3HTHE  ARRAY  OID  NOT  CONVERGE  IN  THE  SPECIFIED,  ) 
*/4X»  1 6HNUM  PER  OF  CYCLES  > 

WRITE!  6,  140  ) 

140  FORMAT!  // 4 X,  1 9HTHE  SMOOTHEO  MATRIX,  //  > 

00  150  T  =  1,  NROW 

WRITE!  6,  70  )  !  G!I,J>,  J  =  1,  NCOL  > 

150  CONTINUE 
STOP 

160  WRITE!  6,  170  ) 

170  FORMAT!  //4X,  31HAT  LEAST  ONE  WEIGHT  IS  NEGATIVE  > 

STOP 

ISO  WRITE!  6,  1 90  > 
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FORMAT (  //RXf  29HR0US  AND/OR  COLUMNS  EXCEED  20  > 

STOP 

END 
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•SUBROUTINE  SMOOTH!  NROW*  NCOL*  X*  W*  NCYCLE*  G*  ICYCLE*  I F  AULT  ) 

SUBROUTINE  TC  ORDER  A  2  DIMENSIONAL  ARRAY  USING 
AN  ALGORITHM  OF  DYKS  TRA  AND  ROBERTSON!  19R2  ). 

THE  PRDERING  IS  DONE  SO  THAT  THE  RE  GRFSS ION 
FUNCTION  IS  INCREASING  IN  EACH  INDEPENDENT  VARIABLE. 

DIMENSION  X!20*20>*  WT!20*20>?  G!20*20>*  C!20*20> 

DIMENSION  R!20*20>*  ROWSI201*  COLS!20>*  WCIGHTl^P) 

DIMENSION  ORDER !20) «  GC HECK ! 2 0 « 20 > *  WI20.20) 

DATA  Ere/  i.OE-Of,/ 

IFAULT  =  0 

CHECK  THAT  ROUS  AND  COLUMNS  DO  NOT  EXCEED  20. 

IF  !  NPCW  • GT «  20  .OR.  NCOL  .GT.  20  )  GOTO  120 

CHECK  THAT  WEIGHTS  ARE  POSITIVE*  OR  ZERO. 

DO  5  I  =  1,  NROW 
DO  5  J  =  1,  NCOL 

IE  !  W!  I*  J  )  ,LT.  0.0  )  GOTO  110 
UT ! 1  * J )  =  W  !  I  *  J  ) 

IE!  WT ! I*  J  >  .GE.  EPS  )  GOTO  5 
UT!I*J)  =  C. 00001 
IFAULT  =  A 
CONT  INLE 

INITIALIZE  R  AND  C  TO  ZERO*  AND  SET  UP  WORKSPACE. 

00  20  I  =  1 ,  NROW 
DO  10  «  :  ],  NC CL 
G!  I*  J  >  =  X!  T*  J  ) 

GCHECK!  I,  J  )  =  X!  I,  J  ) 

C !  I «  J  )  =0.0 
R<  I*  W  )  =  0.0 
CONT  INUE 
CONT INUE 

INITIALIZE  COUNTER  E OR  NUMBER  OE  CYCLES. 

ICOUNT  =  c 

smooth  over  rows. 

UCOUNT  =  0 

DO  b  0  I  =  1  *  NROW 

DO  3  0  v.  =  1  *  NC  CL 

ROWS  !  J  )  =  G!  I,  J  )  -  R<  i,  j  > 

WEIGHT!  J  )  r  WT!  I*  J  ) 

CONT INUr 


CALL  PAV!  NCOL*  ROWS*  1*  WEIGHT,  ORDER  > 


KCOUNT  =  0 

00  AR  U  =  1 «  NCOL 

R!  I*  U  >  =  ORDER!  J  )  -  ROWS!  J  ) 

G<  I*  J  )  =  ORDER!  J  ) 

IF  !  ADS!  G!  I*  J  )  -  GCHECK!  I,  J  )  )  .LT. 
*  KCOliNT  r  KCOUNT  ♦  1 
GCHECK!  I*  U  )  =  G!  I,  J  ) 


EPS  ) 


CONT  INLE 

DETERMINE  IF  thfre  IS  NO  CHANGE  IN  THE  ITH  ROW 
FRO”  THE  PREVIOUS  ITERATION. 

IF  (  KCPUNT  • EQ •  NCOL  )  JCOUNT  =  JCCUNT  ♦  1 
CONT  INU  R 

ICOUNT  -  I C  CUNT  ♦  1 

IF  <  I C  CUNT  .  EQ.  1  )  GOTO  5E 

DETFRM  INF  IF  THERE  HAS  BEEN  NO  CHANGE  IN  ALL 
ROW*'  FROM  THE  PREVIOUS  ITERATION. 

IF  (  JCCUNT  .EO.  NROU  )  GOTO  90 


SMOOTH  OVER  COLUMNS • 

LCOUNT  =  0 

DO  80  J  =  I*  NCOL 

DO  60  I  -  1  .  NROU 

COLS  (  I  >  =  G<  If  J  >  -  C<  It  J  > 

WEIGHT  <  T  )  =  WT<  I,  J  ) 

CONTINUE 

CALL  P A  V(  NROUf  COLSf  If  WEIGHTf  0&CER  > 

MCOUNT  =  n 

DO  70  I  =  If  NR  CW 

C(  It  J  )  :  ORB  E  R (  I  )  -  COLS  C  I  ) 

G<  If  J  )  =  ORDER!  I  ) 

IF  (  APS!  G(  It  J  )  -  GCHECKC  It  J  )  )  .LT.  EPS  > 

*  MCOUNT  =  MCOUNT  ♦  1 
GCHE CK (  It  o  )  =  G<  It  J  ) 

CONTINUE 

DETERMINE  IF  THERE  IS  NO  CHANGE  IN  THE  ITH  COLUMN 
FRO*  THE  PREVIOUS  ITERATION. 

IF  (  MCOUNT  .EQ.  NROW  >  LCOUNT  =  LCOUNT  ♦  1 
CONT  INUE 

ICOUNT  =  I C  CUNT  ♦  1 

DETERMINE  IF  THERE  HAS  BEEN  NO  CHANGE  IN  ALL  COLUMN 
F R 0 v  THE  PREVIOUS  ITERATION. 

IF  (  LCOUNT  •  FQ  •  NCOL  >  GOTO  90 

CHECK  IF  NUMBER  OF  CYCLES  HAS  BEEN  REACHED . 


IF  (  NCYCLE  .EQ.  ICOUNT  )  GOTO  100 
GOTO 

1CYCLE  =  1 C  CUNT 
RETURN 

ICYCLE  =  ICOUNT 

IFAULT  -  1 

RETURN 

IFAULT  =  2 

RETURN 

IFAULT  =  .1 


RETURN 

END 

SUBROUTINE  PAV<  K.  X,  I  ORDER*  U«  FINAL*  ) 

SUBROUTINE  TO  APPLY  POOL  ADJACENT  VIOLATORS  THCOREM . 

DIMENSION  P  X ( 20  )  »  W'T(20>,  PW(20>*  X<20>»  U < 2 0  )  * F I N A  LX ( 2 C ) 
DIMENSION  w  1(20) 

DATA  EPS  /  1.0E-06  / 

SET  UP  UCRKSPACE. 

DO  1  0  T  =  1 ,  K 
r  X  <  I  )  =  X  <  I  ) 

IF<  10»PER  .EG.  0  )  F X <  I )  r  -FX(I) 

UT ( I  )  r  U ( I  ) 

PW(I  )  “  W  T  (  I )  *  FX(I) 

W1  (  I  )  r  U  C 1  ) 

CONTINLE 
IBEL  =  K  -  1 
I  =  0 

I  -  I  ♦  1 

I F (  I  .GT.  IPEL  )  GOTO  50 

II  =  I  ♦  1 

determine  if  pooling  is  required. 

IF<  (  F  X  <  I )  -  FX(I1)  >  .LE.  EPS  )  GOTO  30 

POOL  THL  A  OJ  A  CENT  VALUES. 

PU(I  )  =  PV  <  I)  ♦  PW<  ID 

uki  >  =  w i  ( i >  w i c  ii) 

FX(I  )  -  PW(  I)  /  UKI) 

IP  =  I  •»  1 

I F  (  IP  .GT.  IBEL  )  GOTO  45 
DO  4  0  J  =  IP,  IBEL 
J1  =  J  ♦  1 
PU(J)  =  PU( vl> 

WKJ)  =  VI  l  (  J1 ) 

FX (J)  =  FX( Jl) 

CONT  INUE 
IBEL  =  IBEL  -  1 
GOTO  35 
I  BEL  1  -  IBEL 
ICOUNT  r  0 

IF(  IBF.L1  .LE.  0  )  GOTO  70 

DETERMINE  IF  ALL  VALUES  ARE  ORDERED. 

DO  60  L  =  1«  IPEL1 
LI  =  L  ♦  1 

IF  {  (  FXCL)  -  FX<L1)  )  .LE.  EPS  )  ICOUNT  =  ICOUNT  ♦  1 
CONT INUE 

IFt  ICC  UN  T  .NE.  IBEL1  )  GOTO  20 
RECOVER  FINAL  ORDERED  VALUES. 


FINALXC1)  =  F  X  < 1  ) 


1=2 

IFC  L  *GT.  K  )  f-OTn  loo 
WEIGHT  =  WT  <L  >  ♦  UT(I) 

IF<  WEIGHT  .GT.  W1(J>  )  GOTO  9C 
FTNALXCO  =  F  X  <  J  ) 

WT ( I  )  =  WT (  I)  ♦  WT(L) 

L  =  L  ♦  1 
GOTO  8 C 
J1  =  J  ♦  1 
FINALX(L)  =  FX(Jl) 

I  =  L 
vi  =  J  +  1 

L  =  L  ♦  1 
GOTO  8T 

IF (  I  OR  HE  R  .EO.  1  >  RETURN 
00  110  I  =  l,K 
FINALXU)  =  -FINALX(I) 

CONT  irjur 
RETURN 
ENG 
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